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Abstract 

We investigate stationary nonequilibrium states of systems of particles moving 
according to Hamiltonian dynamics with specified potentials. The systems are 
driven away from equilibrium by Maxwell demon "reflection rules" at the walls. 
These deterministic rules conserve energy but not phase space volume, and the 
resulting global dynamics may or may not be time reversible (or even invertible). 
Using rules designed to simulate moving walls we can obtain a stationary shear 
flow. Assuming that for macroscopic systems this flow satisfies the Navier-Stokes 
equations, we compare the hydrodynamic entropy production with the average rate 
of phase space volume compression. We find that they are equal when the velocity 
distribution of particles incident on the walls is a local Maxwellian. An argument 
for a general equality of this kind, based on the assumption of local thermodynamic 
equilibrium, is given. Molecular dynamic simulations of hard disks in a channel 
produce a steady shear flow with the predicted behavior. 



1 Introduction 

Stationary nonequilibrium states (SNS) of macroscopic systems must be maintained by 
external inputs at their boundaries. Since a complete microscopic description of such 
inputs is generally not feasible it is necessary to represent them by some type of modeling. 
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However, unlike systems in equilibrium, which maintain themselves without external 
inputs and for which one can prove (when not inside a coexistence region of the phase 
diagram) that bulk behavior is independent of the nature of the boundary interactions, 
we do not know how different microscopic modeling of boundary inputs, representing 
fluxes of matter, momentum or energy, affects the resulting SNS. 

What is observed experimentally is that in regimes close to equilibrium, when the 
fluxes are small, the bulk macroscopic behavior is determined by the unique solution of 
the hydrodynamic equations, with specified boundary conditions on the hydrodynamical 
variables such as density, temperature and fluid velocity [1]. The situation may change 
dramatically however as soon as the driving forces become sufficiently large for this 
solution to lose stability. We can then have the formation of coherent structures, such as 
rolls or hexagons, whose pattern is influenced by details of the boundary conditions [2]. 

Even in the absence of hydrodynamic instabilities, e.g. in passive heat conducting 
systems or fluids in regimes of laminar flow, the SNS generally have very long range mi- 
croscopic correlations, with slow power law decay, which can be measured experimentally 
[3]. This raises the possibility that even in regimes of hydrodynamic stability the mod- 
eling of the boundary inputs may have global, albeit subtle, effects on the nature of the 
SNS. In fact it is known that, even in the near equilibrium regime, different modelings of 
the external drives produce very different types of microscopic measures of (what appears 
to be) the same macroscopic SNS. Thus, stochastic drives, such as thermal boundaries 
in which the particles acquire, following a collision with the walls, a specified Maxwellian 
velocity distribution generally lead to stationary measures absolutely continuous with 
respect to Liouville measure [4]. The same is true for systems driven by collisions with 
some simple kinds of Hamiltonian infinite particle reservoirs specified by a given distribu- 
tion prior to the collisions [5] . Deterministic thermostatting schemes, on the other hand, 
yield measures singular with respect to Liouville measure [6-9]. 

It is quite possible, even likely, that these great differences in the structure of the 
microscopic SNS (mSNS) do not have any significant effect on the bulk properties of stable 
macroscopic SNS (MSNS). This is what happens for macroscopic equilibrium systems 
which can be described by a variety of microscopic ensembles, e.g. the canonical or 
micro-canonical [10]. Unlike equilibrium, however, it is far from clear at present how 
to characterize essential global features of mSNS. We do expect however that locally 
such system will be close to an equilibrium state, at least when the inputs are confined 
to the boundaries [11]. Indeed, as we shall discuss further later, this property of local 
thermodynamic equilibrium (LTE) holds the key to understanding SNS of macroscopic 
systems in the hydrodynamic regime. 

Questions regarding the nature of mSNS have recently come to the fore due to the 
combination of computer and analytical investigations of SNS with various deterministic 
thermostatting devices [6-10, 12-18]. The simulations have shown that the stationary 
states produced by these drives behave, at least as far as linear transport coefficients 
and other gross properties of MSNS go, in reasonable accord with known experimental 
and theoretical results. In addition, the simulations have found unexpected interesting 
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microscopic structures in these singular measures, e.g. pairing rules for the Lyapunov 
exponents, a formula for large fluctuations in the phase space volume contraction rate, 
etc. 

The rigorous mathematical analysis of such systems has confirmed some aspects of 
the simulation results [9, 14]. This has led Gallavotti and Cohen [15], see also [10, 
16], to postulate what they call the "chaotic hypothesis", based on Ruelle's principle for 
turbulence: U A many particle system in a stationary state can be regarded, for the purpose 
of computing macroscopic properties, as a smooth dynamical system with a transitive 
Axiom A global attractor. In the reversible case it can be regarded, for the same purposes, 
as a smooth transitive Anosov system." 

This hypothesis was shown by Gallavotti and Cohen to imply, for SNS produced 
by reversible thermostatted dynamics, a formula for the fluctuation in the phase space 
volume contraction rate in agreement with the computer simulations [13] and to be 
generally consistent with known results, at least when the driving is not too large [16], 
see also [17, 18]. It also implies, or even presupposes, a strong form of "equivalence of 
ensembles" , for SNS with specified macroscopic flows. This is, as already noted, certainly 
in accord with experience on SNS close to equilibrium where it can be understood as an 
expression of the existence of LTE. It may however also be true more generally, at least 
in some form. 

In the present paper we investigate a new class of models in which the microscopic 
dynamics in the bulk of the system are Hamiltonian and reflection at the boundaries 
are deterministic and energy conserving. This permits us to define a phase space flow 
X = J-'(X), X a point in (a fixed energy surface of) the system's phase space. In this 
respect our model is similar to the bulk thermostatting schemes mentioned earlier [6- 
8]. Unlike those schemes, however, which modify the equations of motion in the bulk 
of the fluid, something which is computationally useful but has no counterpart in real 
physical systems, our model is fully realistic away from the boundaries. In this respect it 
is similar to models in which the driving force is a "boundary layer" of reservoir particles 
or is given by stochastic thermal boundaries [4, 5, 19]. Our combination of realistic bulk 
dynamics and deterministic boundary drives offers a simple model of an mSNS which can 
be investigated by means of classical dynamical systems theory. It will hopefully lead to 
a better understanding of SNS representing real systems. 

Our main conclusion can be summarized as follows: Deterministic boundary driven 
models, reversible or not, accurately represent the bulk behavior of MSNS. There is an 
equality in these models between phase-space volume contraction, and hydrodynamic 
bulk entropy production for SNS of macroscopic systems in LTE. Plausible arguments of 
why this should be so and of how to connect entropy production in nonequilibrium states 
generated by different modelings of the inputs are given in section 8. This is preceded 
by a detailed description of analytic and computer simulation results for specific models 
producing shear flow. A preliminary account of this work is presented in [20]. 
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2 Description of Models 



To make the analysis as concrete as possible we shall consider here SNS representing 
shear flow in a two dimensional system; we imagine this to be the surface of a cylinder 
of height and perimeter L, or a square box with periodic boundary conditions along the 
x axis, i.e. we identify the left and right sides at x — ±L/2. On the top and bottom 
sides of the box, y = ±L/2, are rigid walls at which stand watchful Maxwell demons who 
make the particles reflect according to rules satisfying the following two conditions: 

(i) the reflected velocity is determined by the incoming one and is energy conserving; 

(ii) the particles at the top wall are driven to the right, and those at the bottom wall 
are symmetrically driven to the left. The purpose of rule (ii) is to imitate moving walls 
and thus produce a shear flow in the bulk of the system. The use of a two dimensional 
system and of symmetric rules for top and bottom is for simplicity only and the reader 
is free to imagine instead a three dimensional channel with different reflection rules at 
top and bottom. 

Since reflections at the walls preserve the particle speed, the reflection rule can be 
defined in terms of angles. For a particle colliding with the top wall, let ip and —ip be 
the angles which the incoming and outgoing velocities make with the positive x-axis, the 
direction of the "wall velocity", so that < <p, ip < ir. At the bottom wall, the angles 
are measured between the velocity vectors and the negative x-axis. Then, any reflection 
rule is completely specified by a function ip = f(ip,v), where v stands for the speed of 
the particle. 

For simplicity, we only study here functions independent of v, which are the same for 
both walls so that 



In particular, the identity function f(ip) = ip corresponds to elastic reflections, and 
ip = 7r — ip gives complete velocity reversal reflections. 
One particularly simple reflection rule is given by 



which, for ir/2 < <p < ix clearly satisfies (ii). Under this rule, particles moving 'in the 
direction of the wall velocity', with p> < ip Q , reflect elastically at the wall, while those 
moving 'opposite the wall velocity', ip > <po will reflect straight back, with the velocity 
vector reversed. This rule is non invertible, and will be discussed further in [18], see also 
[23] . Here we shall focus our attention on two invertible rules that we found particularly 
interesting and which we used in our molecular dynamics simulations. They are 



^ = /(¥>), 



< If, 1p < 7T 



(2.1) 




(pxpo 



(2.2) 



ij = (tt + b) - [(tt + b) 2 - V { V + 2b)f 2 



(2.3) 



with b > and 



tjj = op 



(2.4) 
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with < c < 1. We call these the b-model and the c-model, respectively. 

The graph of the function (2.3) is a circular arc terminating at the points (0, 0) and 
(-7T, 7i") lying below the diagonal ip — (p. This function has the symmetry, (p = = 
7T— /(7T— ip), which makes the dynamics time-reversible. Time reversal symmetry which is 
present in Hamiltonian dynamics as well as in the usual Gaussian thermostatted models 
[6-9], means that the system retraces its trajectory backwards in time following a velocity 
reversal of all the particles. This symmetry plays an essential role in some of the analysis 
in [15]. The c-model is not time- reversible. Similarities and differences between the b 
and c models are therefore of particular interest in determining the range of universality 
present in SNS. We shall in fact see that the b and c models behave in a similar way. 

3 The Stationary State 

Let us consider now the time evolution and the nature of the stationary states we might 
expect with our b or c dynamics for a system of N particles with total energy E; average 
particle and energy densities n = N/L 2 , e — E/L 2 . For the sake of concreteness imagine 
the particles to be hard disks with unit mass and unit diameter so E is just their kinetic 
energy. When b = oo or c = 1 the boundary conditions correspond to elastic reflections 
and so, for N greater than one, but less than some jamming value, we expect that, starting 
with any initial measure absolutely continuous with respect to the Liouville measure on 
the H = E surface, the ensemble density will approach (weakly), as t — * oo, the uniform 
density on this surface, i.e. we expect the micro-canonical ensemble of a system of hard 
disks with mixed periodic and reflecting boundaries to be ergodic and mixing. (There is 
actually, in addition to the energy, also a conserved total x-momentum which we fix to 
be zero and ignore). 

When 6 ^ oo or c ^ 1 the rules will clearly produce a drift to the right near the 
top wall and to the left near the bottom wall. This drift should produce, for L large 
compared to the mean free path (7m) -1 , an mSNS representing a system with a shear 
flow [19, 20]. On the microscopic level we expect now that any initial ensemble density 
absolutely continuous with respect to the microcanonical ensemble will evolve, as t — > 00, 
to a stationary measure ft whose Hausdorff dimension^] is less than the dimension of the 
energy surface [6-10]. Such behavior is proven for a single particle, subjected to an 
external force, moving among a fixed periodic array of scatterers [9]. 

We note that when both walls "move" in the same direction, the x-component of the 
total momentum of the system is a monotone non-decreasing function of the number of 
collisions with the walls. Since this is bounded above, the initial ensemble density must 
converge (at least in some weak sense) to a measure whose support is on configurations in 

1 Here the Hausdorff dimension of the measure fi is, in accord with Young [21] and Pesin [22], the min- 
imum of the Hausdorff dimension of subsets of full measure. (Other definitions of Hausdorff dimensions 
of measures are sometimes used, which makes this notion confusing.) The Hausdorff dimension in the 
Young-Pesin version coincides with the information dimension for systems with nonvanishing Lyapunov 
exponents [22]. Note that the support of the measure p, may or may not be the entire phase space. 
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which the particles all move parallel to the x-axis. This situation is rather pathological. 
We expect our system with top and bottom walls moving in opposite directions, to reach 
and stay in an LTE state at least when the shear is not too large. This is consistent 
with fx being singular: even with its Hausdorff dimension of being only a fraction of the 
dimension of the energy surface [10]. 

Assuming that our system will indeed go, for N large, n and e fixed, to an MSNS 
representing a fluid in shear flow, as is indeed seen in our computer simulations to be 
described later, we consider now briefly the purely hydrodynamical description of such 
an MSNS. This is given by the stationary solution of the compressible Navier-Stokes 
equations for the fluid velocity in the x-direction u{y), the temperature T(y) and density 
n(y), in a uniform channel of width L in which top and bottom walls move with velocities 
±u b in the x-direction, have the same temperature T b , and we impose a no slip boundary 
condition [1] 

These equations, which are derived on the assumption of LTE [11], have the form [1, 

6] 

-p(n,T) = 0, 

(3.1) 




= 0. 

Here p(n, T) is the (local) equilibrium pressure of the system at constant density n(y) and 
temperature T(y), 77(71, T) is the viscosity and n(n, T) is the heat conductivity. Equations 
(3.1) are to be solved subject to the boundary conditions u(±L/2) = ±u&, T(±L/2) = 
Tb and fixed average particle density L~ x J L 2 n (v) dy = n. This gives 

p(n,T)=p(n ,T ), du/dy = U/n, - ndT/dy = J{y) = Uu(y), (3.2) 

where uq = n(0), To = T(0), II is the constant x-momentum flux in the negative y- 
direction, J(y) is the heat flux in the positive ^/-direction, and we have used the symmetry 
of the flow about y — 0. These equations can be solved once p, k and 77 are given as 
functions of n and T. The solution will be unique when the average shear, 7 = 2w fe /L, is 
small, see [24]. 

Eq. (3.2) can be integrated further to give 



u(y) = 

so that 



y Jo 77 



dT 1 fr]\ d 2 



, u (3.3) 

dy 2 \kJ dy 



n = ?77, 1 /^=t/ £ -7%y- (3-4) 
L Jo 77(71, T) 
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For dilute gases, tj/k is a constant independent of n and T in which case 

T{y)=T Q - l -(^u\y). (3.5) 

In some cases the y- variation in r\ and k are so small across the channel that we have 
an essentially linear flow regime with, 

u(y)=iy, T(y) = T - (r?/2/t)7V, and IT = 777. (3.6) 

The constancy of p together with the specified average density then determines n(y). 

4 Entropy Production in SNS 

Entropy plays a central role in determining the time evolution and (final) equilibrium 
states of isolated macroscopic systems. Its microscopic interpretation as the log of the 
phase space volume of all micro-states consistent with a specified macroscopic descrip- 
tion, was well understood by Boltzmann and the other "founding fathers" of statistical 
mechanics, although there is still much fuzziness and outright confusion surrounding the 
subject. We refer the reader to [25] and references there for a discussion. 

The role of entropy and/or entropy production in SNS is also very important although 
much less clear. By their nature truly SNS cannot occur in an isolated finite system evolv- 
ing under Hamiltonian (or quantum) dynamics — the only truly stationary macroscopic 
state for such a system being the equilibrium one. The situation can be different for 
ah initio infinite systems [26], but we shall not discuss that here. We shall instead de- 
scribe now various aspects of entropy production in the simple SNS corresponding to 
stable shear flow in finite systems considered here. We will then discuss in section 8 the 
connection between them and what they teach us. 

a) Hydrodynamics 

The hydrodynamic entropy production a per unit volume in our stationary system is 
given by the Onsager form [17]; see in particular chapter 14 in Balian [1] and section 2.2 
in ref. [7] 




where we used (3.2) in the second equality. The total hydrodynamic entropy production 
R due to the dissipative fluxes in the steady state is then 

R = I adr= [ [Uu/T]ds= f j b /Tds 

J Volume J Surface J Surface 

= J b /T b = 2L 2 B /T b = L 2 B 7 /T b (4.2) 

where j b is the heat flux per unit length and J b the total flux to the walls and we have 
taken the channel to be of length L with periodic boundary conditions in the x-direction. 
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Eq. (4.2) is interpreted in the macroscopic formulation of irreversible thermodynamics 
[1] as an equality, in the stationary state, between the hydrodynamic entropy produced 
in the interior and the entropy carried by the entropy flux, equal to Jb/Tb, to the walls 
of the container. To maintain such a steady state in an experimental situation requires 
external forces acting on the walls to make them move with velocities ±Ub- The work 
done by these forces, |Hu b | per unit wall area, is converted to heat in the bulk of the 
fluid by the viscosity and then absorbed by the walls acting as infinite thermal reservoirs. 
The steady state hydrodynamic entropy production in the system, R, is also carried to 
the walls by this heat flux. If we imagine the walls as "equilibrium" thermal bath at 
temperature T b then R is equal to the rate of their entropy increase (dS cq = dll /T): note 
that we are assuming here that there is no slip between the temperature of the fluid at 
the walls and the temperature of the walls, 
b) Microscopic: Stochastic Reservoirs 

The existence of macroscopic steady states, satisfying the compressible Navier-Stokes 
equations (3.1), can be proven in suitable scaling limits, by starting from the Boltzmann 
equation [24]. In such analysis the walls are typically modeled by stochastic thermal 
boundaries; following collisions with the walls particles have a Maxwellian velocity dis- 
tribution with mean ±Ub and temperature TJ,. Going beyond the mesoscopic description 
given by the Boltzmann equation it is expected that such thermal boundaries will pro- 
duce similar SNS for general fluid systems which will be described, on the microscopic 
level, by a stationary measure on the phase space having a density, p, absolutely con- 
tinuous with respect to Liouville measure [5]. In fact it is possible to show, for systems 
in contact with a thermal reservoir at temperatures T, that the total "ensemble entropy 
production" 

K(t) = S G (t) + (J)/T (4.3) 

is non-negative [5]. Here, 

S G (t) = S G (fi(X,t)) ee J n(X,t)\ogn{X,t)dX, (4.4) 

is the system's Gibbs entropy, dX is the Liouville volume element in the phase space, 
and (J) is the ensemble average of a phase space function J(X) representing energy flux 
to the reservoir. At the same time the rate of change of the mean energy in the system 
is given by 

j f J tiHdx = (H) = —(J) + (W). (4.5) 

where (W) is the average mechanical work done on the system by some external force, 
e.g. one produced by moving rough walls of the system in a channel, c.f. [19]. 

In the stationary state obtained in the limit t — > oo, /i = p, so (H) and S G are 
constant with (W) = (J) and TZ = TZ = (J)/T = (W)/T. Hence if we identify (J) with 
Jb then TZ is equal to R given in (4.2). For a system in contact with only a thermal 
reservoir, the stationary state is the equilibrium one and TZ(t) = d/dt[S G — (H)/T] — > 
as t — > oo, see [27] and section 8. 
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c) Microscopic: Deterministically Driven Systems 

Let us turn now to our models where the flow is deterministic and collisions with the 
boundaries conserve energy. It is not clear at all a priori what should now correspond at 
the microscopic level to the hydrodynamic entropy production in our system. Following 
the work in refs. [7-9] , we note that for a deterministic flow in the phase space the rate 
of change of the systems Gibbs entropy defined in (4.4), is given by 

S G (t) = J n(X,t)(divX)dX = -M(t), (4.6) 

This vanishes for an isolated system evolving according to Hamiltonian dynamics for 
which divX = 0, but not for dynamics which does not conserve Liouville volume. Fur- 
thermore, since we expect the stationary measure for our dynamics, /t, to be singular 
with respect to Liouville measure, we will have Salt) — > — oo. At the same time since 
the convergence of n to jl is in the weak sense, we might have a non vanishing limit 

S G {t) -> j /i(X)(divX) dX = —M (4.7) 

We can then interpret M, the average compression rate of phase space volume per unit 
time as the "measure entropy production" in the stationary state. The existence and 
negativity of the limit in (4.7) was proven in [9] for a simple bulk thermostatted model. 
The non-negativity of M for a stationary ft is proven in a suitable general setting by 
Ruelle [23]. 

The behavior of Sa(t) in the driven deterministic case is to be contrasted with the 
case of stochastic reservoirs considered earlier where the stationary measure has a smooth 
density with Sa(t) — > as t — > oo, while TZ(t) — > 71 = R, the positive hydrodynamic 
entropy production in the MSNS. Now in the bulk thermostatted models, [7-9], the 
equations of motions are such that M is automatically equal to the ensemble average 
of microscopic quantities which can be identified with thermodynamic forces and fluxes 
such as appear in the macroscopic entropy production. This is not the case for the models 
considered here. We have no a priori prescription of u or T anywhere in the system and 
phase space volume gets compressed only at collisions of a particle with the wall: the bulk 
dynamics being Hamiltonian. We therefore need to investigate here the relationship, if 
any, between M and R for our system. Unfortunately, a direct computation, using only 
the given dynamics, is totally out of reach of our present mathematical abilities. What 
we shall do instead in the next section is to make some reasonable assumptions on the 
nature of the microscopic SNS in the limit when our system becomes of macroscopic 
size. It will turn out that these assumptions, which are satisfied for a system in LTE, 
lead to an equality between M and R. This will be checked and confirmed by computer 
simulations in section 7. We will then argue, in section 8, that such an equality holds in 
general when the macroscopic system is in a state of LTE. 

Remark. It might be feasible to carry out such a rigorous analysis of our model within 
the context of the Boltzmann equation, in analogy to what is done for stochastic walls in 
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[24] . Numerical simulations on models b and c using the Direct Simulation Monte Carlo 
method for simulating the Boltzmann collision term inside the channel are now being 
carried out [26]. The results appear consistent with those in section 7. 



5 Calculation of M in the Hydrodynamic Regime 

To obtain the rate of compression M for our system, we observe that our dynamics 
preserves phase space volume except at collisions of a particle with a wall. Since these 
collisions take place "instantaneously" we can compute the compression occuring at a 
single collision — ignoring the rest of the particles. The compression is then just equal 
to the ratio of the "outgoing" one particle phase space volume (dx' dy' dv' x dv' y ) to the 
incoming one (dx dy dv x dv y ) in a time interval dt containing the collision. A little thought 
shows that dx' = dx and \dv' x dv' y /dv x dv y \ = \v' dv' dip\/\v dv d(f>\ = \dip/d(p\ where and 
ip are the incoming and outgoing angles, related by ip = f (</>). Similarly, \dy'/dy\ = 
\v'ydt/v y dt\ = | sinip/ sinc/>|. Hence in every collision between a particle and the wall the 
phase space volume is changed by a factor 



sin ijj dip 



sin p dip 



sin<£> 



d cos f(ip) 



dcosp 



(5-1) 



where / defines the reflection rule (2.1). The phase volume will be reduced or increased 
depending on whether (5.1) is smaller or larger than unity. 

The mean exponential rate of compression of the phase volume per unit time is then 
given by 

M = -2(N c \og[f'(p)smf(p)/smp]) fl (5.2) 

where N c (tp) is the flux of particles entering a collision with the top wall at angle ip. The 
factor 2 comes from summing over top and bottom walls and the average is taken with 
respect to the stationary measure fi. 

An exact evaluation of M given in (5.2) is currently far beyond our abilities. To pro- 
ceed further we assume now that in the hydrodynamic regime corresponding to L ^> /, 
and 7/ <C 1, where I ~ (7m) ~ x is the mean free path between particle-particle collisions, 
the density p(vi,v 2 ), of the velocity vectors, v = (vi,v 2 ), of particles entering a colli- 
sion with the walls, is Maxwellian with the (to be determined) mean value (v, 0) and 
temperature T w for the top wall, ((—v, 0) and T w for the bottom wall). That is 

p( Vl ,v 2 ) = (2nT*y^v 2 e W (- K ~5 + ^ 1 , > (5-3) 



2T, 



w 



It will be convenient to rewrite the density (5.3) in polar coordinates, V\ — r cos 9 and 
v 2 = r sin 9, 
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and denote by (F) the average of any function F(r,9) with respect to (5.4). 
The average momentum transfer 'from wall to wall', II, is now given by 

fl = n c (Avi) = n c (v liOUt - v i;m ) (5.5) 

where n c is the 'collision rate', i.e. the average number of collisions with the top wall per 
unit time per unit length, which will have to be computed later. In our polar coordinates 
(r, 9) this value is 

ft = n c (r(cos f{6) — cos 6)) 

= c -^r / / (cos/(0) -cos#)r 3 sin# 

( v 2 + r 2 — 2vr cos 6\ in 1 . . 

x exp I — I d9dr (5.6) 

In order to keep the "velocity" of our walls from growing with L when L becomes 
macroscopic, which would certainly take the system away from local equilibrium, we need 
to consider situations in which ft, like the hydrodynamic II in sec. 3 is of order (yj). This 
requires that the reflection rules (2.1) be close to the identity, i.e. we put 

i, = f(tp) = up + Shiip) + 5 2 f 2 {uj) + o(5 2 ) (5.7) 

with fi, fi some fixed functions on [0, n]. In fact, as seen from (5.6), 5 has to be of order 
0(L _1 ). For our b-model (2.3) we can set 5 = 1/6, and then (5.7) will have the form: 

ijj = ip- Sip(ir -<p) + 5 2 ip(n - ip) 2 + o(5 2 ) (5.8) 

For our c-model (2.4) we can set 5 = 1 — c, and so 

if) — (p — Sip (5.9) 

(Our calculations however are not restricted to these two models. ) 
Expanding ft in 5 gives 



IT = n c {r(-sm6-5f 1 (6)+0(5 2 )) / 

n c 5 f°° r w 3 . 2 „ ( v 2 + r 2 — 2vr cos^" 



f°° r t tn\ 3 • 2/1 \ V + r ZvrC0SU \ jnj 

Jo Jo ^ ^ XP I 2T J 



2nTl' 2 Jo Jo \ 
+ 0(5 2 ) (5.10) 

The last double integral depends on the so far unknown parameters v and T w , and we 
denote it by I(v,T w ). 

We next assume that the steady state will indeed correspond to a shear flow described 
(on the average) by the hydrodynamic equations in sec. 3, and identify the II with II and 
T w with Tb there. Since we are however not given T& and uj, they have to be determined 
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from the given data, assuming the system to be in LTE. Thus the determination of T w 
is now done on the basis that the system has a fixed energy E so that 

j £n(y)[T(y) + \u 2 (y)]dy = E/L 2 = e, (5.11a) 

while the mean horizontal velocity of the particles near the wall can be taken as the 
mean between ingoing and outgoing velocities, 

U b = (ui )OU t + Ul,in) /2 = v + 0(5) (5.116) 

The solution of (3.1) which determines u(y), T(y) and n(y) in terms of u^, 7& and n will 
now be determined entirely by the a priori given e, n and the rule /(<£>), via (5.10) and 
(5.11). The computation becomes straightforward in the linear approximation (3.6), see 
eqs. (6.6)-(6.9) in the next section. 

We study now the case 5^0 and L = a/5 with some fixed a > 0. Then, the wall 
velocity u b is proportional to II L and thus does not vanish as 5 — > 0. Using (4.2) and 
(5.10) the hydrodynamic entropy production, R, is then, in this limit, given by 

R=- 2 ™ C l 2 mT w ) + o(l) (5.12) 



To obtain the compression rate M in this limit we expand (5.2) in 5 using our ansatz 
(5.3). This gives, upon replacing N c by 2n c L, 

M/(2n c L) = -(log[/V) sin /(</?)/ sin </?]> 

= - (AH cos ip I sin <p) 5 -(f[(<p)) 5 + o(5) (5.13) 

Now, integration by parts yields 

POO f"K 

(fM) = - / fi(e)pe(r,9)d9dr 
Jo Jo 

where pe stands for the partial derivative of the density (5.4) with respect to 9: 

. n . cos9 vrsin9 
Pe(r, 9) = ~^p(r, 9) - — -p r, 9) (5.14) 
sin v l w 

Combining (5.13) and (5.14) then gives 

M = -25Ln c vT^ (f^r sin tp) + o(L5) = - c I(v,T w ) + o(l) (5.15) 

V 27T T w 



A shorter way to get (5.15) is to rewrite the middle in (5.13) as 

log (l + MO) )=-s( l^)^-) + o(6) (5.16) 
a cos ip J \ a cos ip ' 
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and then do an integration by parts using cos6> as a variable. The leading term in 
(5.16) will be recognized as corresponding to (4.6) for the continuous time action of the 
thermostats, see Appendix. 

The leading term in the expansion of M is thus exactly the same as in that of R, 
hence M and R become equal in the hydrodynamical limit, L — > oo. The essential re- 
quirement for the equality is the validity of (5.14). (This is weaker than (5.3), permitting 
multiplication of the Maxwellian there by an arbitrary function of r, but we do not know 
of any physically reasonable non Maxwellian p which would satisfy (5.14).) The equality 
between R and M in the hydrodynamical regime is thus a nontrivial consequence of our 
(local equilibrium) assumption (5.3). As already noted this is an important difference 
between our models and the bulk thermostatted models of [7, 8, 9]. The relation between 
phase-space volume compression and what looks like entropy production is so built into 
the structure of the dynamics of the latter that there equality holds, essentially by def- 
inition, even for systems consisting of just of one or a few particles which are certainly 
not in local equilibrium. This is not the case here. The equality fails for systems with 
too few particles to be well described by hydrodynamics, as is seen in the next section. 



Before presenting our numerical simulations, which were obviously done at finite (and 
not so large) L, we consider the consequences which can be drawn from our assumption 
(5.3) when 5^0 and L is (relatively) large but held fixed. Now the generated shear 
flow in our system will only be approximately described by the hydrodynamics equation 
(3.1). Furthermore u w as well as R and M vanish as 5 — > 0. Interestingly enough there is 
now a difference between the b and c models with M/R remaining finite for the b-model 
as 5 — > (close to unity for moderately large L) , while becoming infinite for the c- model. 
As we shall see, however, this is connected with the details of the c-model, rather than 
with its lack of time reversibility. 

To compute R and M in this case it will be necessary to expand various quantities 
up to second order in S. Thus, 



This expansion is of course meaningful only if the coefficients of 5 and 5 2 are finite which 
requires that /i(0) = fi(n) = and that f\ has finite derivatives at both and n. Our b- 
model satisfies these assumption, because fi(<p) = —p(7T — p), see (5.8). The computation 
of M for the c-model which does not satisfy f\ (ir) = will be done separately. 



6 Small S regime 



M/(2n c L) 



+ 



(log[/'M sin /(</?)/ sin <p]) 
(fi(<p) cos <p/ simp) 5+ (f[(p)}5 
(f 2 (<p)co8<p/sm<p) 6 2 + (ti(p>))6 2 



(6.1) 
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In order to compute the coefficients in the expansion (6.1) we need to expand the 
density (5.4) which using (3.2) or (3.6) depends on 5 through the quantities v = 0(5), 
T w = l/2 + 0(5 2 ) to first order in 5. This gives 

p(r, 9) = 2ir- 1/2 r 2 e- r2 sin 9 + A-K~ 1/2 vr z e~ r2 sin 9 cos 9 + o(5), (6.2) 

M = (d + C 2 )v 2 + o(5 2 ) (6.3) 

with 

c '=*'( i -^r (6 - 4) 

C 2 = 2 -^j- (l - -W ( I" /i(«) S m 2 « deY 2 - sin" 1 0+(/«fl)) J staff] <W (6-5) 

n c L \ n c LJ \Jo J Jo 

We assume now that for small 5 and moderate L, of the kind used in our simulation, 

the system is reasonably well described by a linear shear flow, Eq. (3.6). We can then 

find the relation between v and 5, 

Lnr6 f ' V ~Y £ M9) sin 2 9 d9 + 0(5 2 ) (6.6) 



2t7a/7t V n c L, 



Similarly we get 



n = r r h{9)r 3 sin 2 9- e^ 2 d9dr + 0(5 2 ) 

-J ix Jo Jo 

" r() r h{9) sin 2 9 d9 + 0(5 2 ) (6.7) 
Jo 



^/n Jo 

R = ^ = uW( 1 _j,y\ (f maaHdS y +0(P) (6 .s) 



and 

R = 2 - 

T w nr] \ n c LJ \Jo 

The expression for M in (6.3) in terms of 5 and L is 



TTT] 

1 t ,2 V 



+ -Ln c S 



n,L 



r\fl(9) sin" 1 9 + (f[(9)) 2 sin 9} d9 + o(5 2 ) (6.9) 
Jo 



Comparing (6.8) to (6.9) we see that the term of M proportional to L 2 coincides with 
the corresponding term of R while the terms of order L differ. The difference between 
M and R thus gets relatively small as L increases with the value of M remaining slightly 
larger than the value of R as 5 — > 0. This is in agreement with the computer simulations. 

We now calculate M for our c-model where, according to (5.9), fi((p) = —<f and 
/ 2 (<^) = 0. The expansion (6.1) now has a singular term and so we need to use the 
original formula (5.2) 



M = -2Ln c (ln(l -5) + (lnsin((l - 6)ip)) - (lnsin^) 

= -2Ln c (I + h-I 2 ) (6.10) 
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where Jo = — 5 + 0(5 2 ) and the two other terms involve integration with respect to the 
density (6.2). It will be sufficient to integrate only with the zeroth order term there 
since corrections will, as before, produce terms of order 5 2 , which we can now disregard 
compared to the leading term. We find 

I 2 = (lnsiny?) = - f sm9 \nsmddd = ln2 - 1 (6.11) 
2 Jo 

and 

h = (lnsin((l -6)<p)) = I r sm9\nsm((l- 5)0) d9 (6.12) 

2 Jo 



This yields: 



7T 2 

h = ln2-l + 5+ — 5 2 \n5 + 0(5 2 ) (6.13) 



Combining the expansions for I , Ii and J 2 gives 

M = -2- l 7c 2 Ln c 5 2 \n5 + 0(5 2 ) (6.14) 
The corresponding expression for R has no singular terms 



7i 3 L 2 n 2 



R= " ~ '" c (1- -\)s 2 + o(5 2 ) (6.15) 



8r] V n c L 
This gives 

M/R = const • ln(l/5) + 0(1) -> oo (6.16) 

In other words, for our c- model, when L is finite and 5^0, the quantities M and R 
have different rates of convergence to zero! 

Remark. We also calculated the second order term in the expansion of M: 

M = -2- 1 vr 2 Ln c 5 2 In 5 + C 2 5 2 + 0{5 3 ) (6.17) 

where 

C 2 = — ^ 11.9202 n c L (6.18) 

7 Numerical results 

Computer simulations of our b and c models were performed with iV = 100 and iV = 200 
hard disks of unit mass and diameter a — 1. We kept the volume fraction occupied by 
the disks, nn/A, equal to 0.1. Our system was thus in the dilute gas phase. For N = 100 
this corresponds to L — 28.0 and for iV = 200 to L — 39.6. The mean free path / at this 
density is about 2.3 [27] so L/l ~ 17 for the latter. (Discrepancies between results of the 
simulations, and of hydrodynamics, can therefore be expected, a priori, to be of order 

{2u b /n).) 
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The positions of the particle centers were chosen randomly (without overlap) in a 
domain \x\ < L/2, \y\ < (L — l)/2. The initial velocities were generated randomly with 
subsequent normalization of the total kinetic energy so that 2E = £ v? = N. 

There is an obvious instability in the dynamics, and the round-off errors accumulate 
at an exponential rate. In fact, the system loses its memory completely after every 100- 
200 collisions in the box. The dynamical meaning of computer simulations of unstable 
dynamical systems like ours is a difficult question, which has been discussed recently in 
the literature [28]. One way to think of round-off errors is as small perturbations of the 
deterministic trajectory of the phase point made at every collision. So, instead of a true 
trajectory we can only track a perturbed one, or a pseudo-trajectory. Then the question 
is - what are the quantities measured by averaging along such pseudo-trajectories? 

One possible answer is given by the so-called shadowing lemma [29]. It says that 
for smooth hyperbolic systems every pseudo-trajectory is shadowed (approximated) by a 
real trajectory. The distance between the two trajectories is of the order of the computer 
accuracy. This may justify averaging over computer-generated pseudo-trajectories for 
smooth hyperbolic systems although it certainly does not guarantee that these finite 
time averages represent typical behavior [28]. Moreover, our dynamics are not smooth, 
and hyperbolicity cannot be easily established. Trust in the results of our simulations, 
like all others done on such unstable dynamical systems therefore relies mainly on faith 
in some kind of typicality resulting from the (pseudo) random effects of the roundoffs 
[28]. 

In any case to prevent the system from leaving the energy surface due to round-off 
errors, the total kinetic energy was renormalized after every N collisions in the box. 

With the above reservations we believe that the results to be described are statistically 
reliable within one percent. This is based on various checks comparing different runs and 
different levels of computer precision. For each value of b and c we averaged over about 
25,000 collisions per particle with other particles and about 1,200 with the walls. We also 
changed the computer precision from single (7 accurate decimals) to double (14 decimals) 
and the so called long double (19 decimals) to make sure the results were stable. Most 
of the programming was done in the C language on an 80486/DX-66 PC and on a SUN 
SPARC station 1000 at the University of Alabama at Birmingham. 

In each run the vertical height, occupied by the centers, L — 1, was divided into 
twenty equally spaced horizontal layers and time averages of the density n(y), mean x- 
velocity u(y) and variances ((v x — u) 2 ), (v 2 ) were taken. We also recorded time averages 
of x-momentum transfer from the walls, II, and the compression rate M. 

Figures 1 and 2 show typical velocity and temperature profiles u(y) for the b- and 
c-models with N = 200 particles. The velocity profiles are almost linear and the tem- 
perature profiles T(y) = \{(y x — u(y)) 2 + v 2 ) almost quadratic, away from the walls, 
consistent with the approximation (3.6). The deviations from linearity near the walls is 
due both to the dependence of rj and koiiT and n which means that we should use (3.3) 
rather use (3.6). Since 77(71, T) ~ y/T this will indeed lead to increases in the slope near 
the wall. There are of course also effects due to the deviations from the hydrodynamic 
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Figure 1: Velocity profiles u(y) for the models b = 10 and c = 0.9. 

limit of 0(2ubl/L), but these are harder to compute. 

We estimated the shear rate du(y)/dy, by the least square fit of a line u(y) = 'jy 
to the experimental velocity profile and used the experimental momentum transfer from 
wall to wall, II, to find the average shear viscosity as defined in (3.4), fj = IT/7. To 
mitigate the problem arising from the nonlinearity of the profile near the wall we used 
only the data in the bulk of the system, from level 4 to level 17, i.e. we excluded the top 
three and the bottom three levels. Table 1 presents the experimental values of the shear 
viscosity computed for the b- and c-dynamics with N = 200 particles. The essentially 
linear dependence of II on S is certainly consistent with (6.6). 



model 


■H-exp 




Enskog f) 


b=10 


6.03 x 10" 


-3 


0.204 


0.209 


b=25 


3.12 x 10- 


-3 


0.213 


0.219 


b=50 


1.66 x 10- 


-3 


0.221 


0.222 


b=100 


0.85 x 10- 


-3 


0.218 


0.222 


b=200 


0.43 x 10- 


-3 


0.216 


0.222 


c=0.90 


4.87 x 10- 


-3 


0.213 


0.214 


c=0.95 


2.82 x 10" 


-3 


0.221 


0.220 


c=0.97 


1.77 x 10- 


-3 


0.206 


0.221 


c=0.99 


0.62 x 10- 


_ 


0.229 


0.222 



Table 1. Computed and theoretical values of 77. 
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Figure 2: Temperature profiles T(y) for the models b = 10 and c = 0.9. 



The shear viscosity of hard disk fluids can be estimated via Enskog's modification of 
the Boltzmann equation [30, 31]. It gives the value 



VE — ''/dilute 



X' 1 + bn + 0.8729x(H 5 



where r/diiute is the value obtained from the Boltzmann equation [30] 



1 / kT 

^dilute = 1.022 • -W — 
Z V vr 



(7.1) 



(7.2) 



for disks of unit mass and diameter. In (7.1) b is the second virial coefficient, b = 7r/2, 
and x is the Enskog scaling factor which is just the equilibrium pair-correlation function 
at contact. We have estimated x by using the pressure equation for hard disk fluids, 



V 



nkT 



71 

1 + 



(7.3) 



and its virial expansion in the number density n, or equivalently the "scaled particle" 
approximation see [33, 34] 



P 



nkT 

In our case bn = 0.2, and we get 



1 + bn + hn 2 + 6 4 n 3 + . . . 



1 - 



7m 



(7.4) 



T]E ~ 1.081 dilute 



(7.5) 
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The value of t]e so computed using the measured mean temperature and density in 
rows 4-17 is shown in the last column of Table 1. The agreement is very good as would 
be expected for the low density system we are considering. 

We next checked the constancy of the pressure p(n,T). To do this it is important 
to note that Eq. (7.3) is valid, even for uniform equilibrium systems, only in the bulk, 
i.e. at distances from the wall large compared to the equilibrium correlation length. For 
our dilute gas system this length is of order a — 1. Near the wall the density becomes 
nonuniform with the density at the wall, n w , equal to p/kT; T the uniform equilibrium 
temperature [30]. Since the hydrodynamic eqs. (3.1) are valid on a scale which is very 
large compared to any microscopic scale this variation in density is not considered there. 
The situation is very different however for our computer simulations where we can expect 
to see these density variations, c.f. [19]. Indeed, the pressure defined as the average y- 
momentum transfer from the top wall in the ^-direction per unit length and unit time 
would be using (5.3), in analogy to (5.6), equal to n c (r(sin/(0) + sin#)). This would 
reduce to n w kT to zero order in 5. 

We present in Table 2 experimental values of n, T, as well as the bulk pressure 
defined in (7.3) in the different layers taking the mean of the values in the layers situated 
symmetrically about the middle. Here the density is given by the average number of 
particles on each level. Since n(y) increases rapidly as we approach the wall, we do not 
have a simple formula for p in the top column. Using n w = p/kT leads to an extrapolated 
density at the wall of about 13.4 and 13 for the b = 10 and c = .9 cases. 







b = 


10 






c = 


0.9 




level 




T 


n 


p 




T 


n 


p 


1 


0.584 


0.399 


11.35 




0.449 


0.439 


10.70 




2 


0.502 


0.420 


10.35 


5.330 


0.383 


0.453 


10.15 


5.722 


3 


0.437 


0.430 


10.13 


5.338 


0.333 


0.459 


10.06 


5.724 


4 


0.374 


0.438 


9.97 


5.337 


0.289 


0.465 


9.94 


5.713 


5 


0.315 


0.443 


9.86 


5.331 


0.244 


0.469 


9.89 


5.712 


6 


0.257 


0.449 


9.77 


5.341 


0.197 


0.471 


9.85 


5.713 


7 


0.198 


0.454 


9.67 


5.341 


0.153 


0.473 


9.84 


5.715 


8 


0.141 


0.455 


9.67 


5.342 


0.111 


0.474 


9.82 


5.720 


9 


0.084 


0.457 


9.63 


5.344 


0.067 


0.475 


9.82 


5.718 


10 


0.028 


0.459 


9.60 


5.340 


0.022 


0.476 


9.80 


5.714 



Table 2. Experimental measurements of the x-velocity, temperature T, density n and 
pressure p on all levels for the b = 10 and c = 0.9 models. 

The agreement between the prediction of the hydrodynamic equation (3.6) and the 
simulation results shows that the Maxwell demon boundary drives indeed set up, in the 
limit L — > oo, an MSNS for shear flow. We therefore computed the values of R and M 
defined for our shear flow by formulas (5.2) and (4.2) with II given by (5.5) experimentally, 
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by using time averages of the appropriate dynamical functions in our simulations. We 
also computed these quantities numerically according to our integral formulas using the 
Maxwellian ansatz (5.4). In both these computations we used the experimental values of 
T w , v w , and n c . The results are presented in Table 3. The last column of Table 3 gives 
the leading term in 5, computed from (5.12), for which M — R. 



model 


M th /M cxp 


Rth/ -Rexp 


M lcad 


b=10 


0.761/0.740 


0.695/0.767 


0.867 


b=25 


0.161/0.157 


0.156/0.162 


0.169 


b=50 


0.0429/0.0417 


0.0422/0.0428 


0.0437 


b=100 


0.0113/0.0110 


0.0111/0.0112 


0.0113 


b=200 


0.00294/0.00285 


0.00289/0.00289 


0.00291 


c=0.9 


0.448/0.444 


0.405/0.432 


0.402 


c=0.95 


0.149/0.148 


0.127/0.131 


0.126 


c=0.97 


0.0641/0.0632 


0.0523/0.0531 


0.0518 


c=0.99 


0.00871/0.00868 


0.00628/0.00632 


0.00625 



Table 3. The experimental and theoretical values of M and R. 



The agreement between the so computed theoretical values of M and R and their 
experimental ones is quite good. This suggests that the integral formulas (5.6) and (5.2) 
with the Maxwellian density (5.4) are quite accurate. 

We also tested directly this hypothesis using a chi-square test and the Kolmogorov- 
Smirnov test, see, e.g., [36]. Both tests accepted the distribution (5.3) for the velocity 
vectors of incoming particle colliding with the wall for various values of b and c. Just to 
determine the sensitivity of our procedure, we tested in the same way the distribution 
of the outgoing velocity vectors (whose Maxwellianity would contradict our assumption 
(5.3)). assumed, and which would have contradicted our assumption (5.3)). As expected, 
both tests frequently rejected this last hypothesis for several values of b and c, indicating 
that the reflection rule (2.1) distorts the velocity distribution in a considerable way even 
for small 5. 
Remarks 

1) Our analysis in section 6 shows that the quantities M and R are of order as 
5^0, with the exception of M for the c-model, which is of order v 2 w \n(l/v w ). Figure 3 
shows the ratios M/v^ and R/v^ as functions of 5. For the b-model they both 'nicely' 
converge to the same positive constant 1.6) as 5 — > 0. For the c-model the ratio 
R/v^j also converges to a number (ss 1.55), while M/v^ apparently grows to infinity, in 
agreement with the prediction in sec. 6. 

2) An interesting question is how the velocity near the wall v w depends on the model 
parameter 5 as the size L is hold fixed and 5 is not very small. Fig. 4 shows the experi- 
mental values of v w versus 5 as 5 varies between zero and .5. A linear regression for small 
5 is very clear for both models. For large S, the function v w (S) increases more and more 
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Figure 3: The ratios M/v 2 w (solid circles) and R/v 2 w (hollow circles) versus 5 = 1/b and 
5 = 1 - c. 
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Figure 4: The velocity v w near the wall versus 5 = 1/b and S = 1 — c. 
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slowly, and apparently has a finite asymptote v w (+oo). This happens for a simple reason 
- the energy balance in the system. Given a linear velocity profile u(y) = 73/ the energy 
balance with the minimal temperature T(y) = gives 7 max = V&/L, so that 

*Vmax = V6/2 « 1.225 (7.6) 

Indeed, the largest velocity near the wall we observed in our simulations (say, with 6 = 
or c = 0.2) was around 1.2. Under these conditions, however, the laminar velocity profile 
breaks down. 
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8 Discussion 



In Sec. 5 we obtained an equality between the phase volume contraction M defined in 
(5.2) for the mSNS produced by our model and the hydrodynamic entropy production 
R defined in (4.2) for the MSNS, in the limit L -> oo, 5^0, with 5L = a fixed. To 
understand the origin of this equality, also found approximately in our computer simu- 
lations, we will analyze here in more detail, the production of entropy in nonequilibrium 
macroscopic systems, discussed in Sec. 4 for SNS. This will make use of formal manipu- 
lations of various expressions for the entropy of such systems whose justification requires, 
at the minimum, the validity of dissipative hydrodynamics, e.g. the Navier-Stokes eqs. 
of section 3, obtained as a scaling limit in going from microscopic to macroscopic de- 
scriptions of our system. It will assume ipso facto the existence of LTE in these systems 
since this is required for the derivation of the hydrodynamic equations [11]. While these 
assumptions are very reasonable for systems not too far from equilibrium, where the 
interactions (collisions) between the particles, tending to bring the system to equilib- 
rium, dominate locally over the external forces and fluxes pushing the system away from 
equilibrium, such results are very far from being proven for systems with Hamiltonian 
dynamics. Even their derivation from the Boltzmann eq. is still incomplete at the present 
time [24]. Given this situation there seems little point in giving any proofs here — even 
for those parts where this may be possible. Instead the analysis should be thought of as 
heuristic and suggestive. 

For a system in LTE with hydrodynamical variables n(r, t), u(r, t), e(r, t), correspond- 
ing to particle density, velocity, and energy density, evolving according to hydrodynamical 
equations, the hydrodynamic entropy, Sh, of the system at any time t is the integral of 
s eq (n, e'), the entropy density in a uniform equilibrium system with densities n and e', 
[1, 11, 25], 

S h = I Scq (n(r),e'(r))rfr (8.1) 

J Volume 

where e'(r) = [e(r) — |n(r)u 2 (r)] is the thermal energy density. For an equilibrium system, 
u has to be independent of r and can therefore be removed by a Gallilean transformation. 

We also have, essentially from the definition of LTE, that macroscopic systems in 
LTE have a local microscopic description which is, to leading order, the same as that for 
a uniform equilibrium system with the same parameters; i.e. if we consider a "small" 
macroscopic volume element around r, its properties will be approximately given by the 
grand-canonical ensemble, assumed to be equivalent to the canonical or micro-canonical 
ensemble, specified by the local values of n, u, e, at time t [11, 34]. The dissipative fluxes 
are then related by the transport coefficients to the gradients of the hydrodynamical 
variables, which are very small on the micro scale. 

Let us call M. the macroscopic state specified by the hydrodynamical variables {n, u, e} 
and let v[X\ M.) be the grand canonical ensemble with local chemical potential and tem- 
perature appropriate to n(r),e'(r). Then we have that Sh is equal to the Gibbs entropy 
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Sq defined in (4.4), corresponding to /i — v, 

S h ~ S G (u) = - J v(X- M) log v(X; M) dX, (8.2) 

see [11, 25, 34]. 

It is furthermore true, as observed by Boltzmann, that the hydrodynamic entropy Sh, 
of a system in the macrostate A4, is equal to the logarithm of the phase space volume, 
T(Ai), associated with the macroscopic state M. = {n, u, e}, see [25], i.e. 

S h ^ S B (M) = log F(M) (8.3) 

(Note that Sb{M) may make sense even when M does not correspond to an LTE state, 
when S h is not well defined. We shall however not consider such cases here, so (8.3) will 
always hold.) 

The equality between the different expressions for the macroscopic entropy of systems 
in LTE given in (8.1), (8.2) and (8.3), depends crucially on the large separation of scales 
between the micro and macro descriptions. It holds to leading order in the ratio of micro 
to macro scales, e.g. l/L defined in sec. 7, and becomes exact, in the sense that it has the 
same limiting value when divided by the number of particles in the system, only when 
the ratio of micro to macro scales goes to zero, the hydrodynamical scaling limit. For 
real systems the equality is only approximate, as are the hydrodynamic eqs., etc. We 
shall however, following Newton in a different but similar context [35], assume that "the 
error will not be sensible; and therefore this ... may be considered as physically exacf 
(italics added) and will therefore treat (8.2) and (8.3) as true equalities. 

Taking the time derivative of (8.1) using the standard equilibrium relations for deriva- 
tives of s eq , we obtain 



s k =l i 

J Volume I 



A dn 1 
T~dt + T 



de d fl 2 



dt dt\2 HU 



dr, (8.4) 



where A(r, t) is the local chemical potential and T(r,t) the local temperature. The 
integrand in (8.4) can be rewritten as = — -j s (r, t) + er(r, t), where j s is the entropy 
flux and o~(r, t) is the local entropy production. The latter can be written in a form 
similar to (4.1) in terms of the full pressure tensor P(r, i), heat flux J(r,t), etc., and is 
non- negative by the second law, see [1, 17, 34]. We thus find, 

S h (t) + I j s (r, t) ■ ds = / <j(r, t) dr = R t (t) > 0, (8.5) 

J Surface J Volume 

For an isolated macroscopic system the integral of j s over the surface vanishes and the 
rate of change of S h (t) is then given by Ri(t). We can thus interpret Ri as the rate at 
which Sb is produced inside the system when it is in LTE. (LTE was implicitly assumed 
in sections 3 and 4a.) Furthermore, since the system is isolated we expect it to approach, 
as t — > oo, a uniform global equilibrium state, with Sh — > •S'cq = (Volume) • s cq (n, e') and 
/.'/ -0. 
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Consider next the hydro dynamic description of an "open" macroscopic system sub- 
jected to external forces and able to exchange heat through its surface with k heat baths 
maintained at temperatures T a , a = 1, . . . , k. (In the shear flow case we could have the 
walls at different temperatures.) We will now have, using (8.5) for the time derivative of 
the system's hydrodynamic entropy, assuming that the part of the system in contact with 
the ath heat bath has a temperature very close to T a , (we are considering for simplicity 
the "no slip" case), 

S h (t)=R l (t)-J2Ut)/T a . (8.6) 

Here J a is the heat flux from the system to the ath heat bath, i.e. J a /T a is the net 
hydrodynamic entropy flux leaving the system through that part of the surface which is 
at temperature T a . Thus J a /T a is analogous to the equilibrium relation dQ = dS/T, 
extended to quasi- stationary processes; see in particular chapter 14 in Balian's book [1]. 
Note that, in accordance with (8.2), the external mechanical forces do not contribute 
directly to the entropy change as they do not, by Liouville's theorem, change the phase 
space volume F(A4) available to the macro state M.. 

When the differences in the T^s and the magnitude of the external forces are not too 
large we expect this system to come to a stationary LTE state in which Sh vanishes, with 

Ri — ► R-i = ^2 Ja/T a \ (8.7) 

in accord with (4.2), for laminar shear flow. 

Let us investigate now the statistical microscopic description of these macroscopic sit- 
uations. Following the usual procedures of statistical mechanics for macroscopic systems, 
we represent their macroscopic state, A4, at time t, by a suitable ensemble density fi(X, t). 
This jj, is assumed to have the property that the hydrodynamic variables M.^ = {n, u, e} 
obtained as expectation values with respect to ^(X, t) of the corresponding phase space 
functions are sharply defined (very little dispersion with respect to fi), vary slowly in 
space and time on the microscopic scales and evolve, on the appropriate macroscopic 
scales, according to the hydrodynamic equations. We will associate to n(X,t) = ^t(X) 
the local equilibrium ensemble density v{X\M. jJit ) and, with some abuse of notation, shall 
set v{X\M. iH ) = v t (X). N.B. While we assume that our system is in LTE, we are not 
assuming here that /j, t — v t . 

We now take the Gibbs entropy Sq of the ensemble density /i t defined in (4.4), and 
split it into two parts, 

S G (fh) = S G (t) = -J lh {X)\og lh {X)dX 

= - J fi t (X)\og[fi t (X)/u t (X)]dX + J fr(X)logv t (X)dX (8.8) 

The first term on the right side of (8.8) is the relative entropy of \x± with respect to u t , 
(which is negative, or zero, by convexity of xlogx). To evaluate the second term we 
use the fact that, by the definition of u, as a locally grand-canonical ensemble density, 
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log v(X) is a function of X whose average is expressible in terms of M. = {n, u, e}. Now, 
by the definition of v t these are the actual hydrodynamical variables in our system so 
that the second term is equal to / vt{X) log i > t(X)dX. Hence the second term in (8.8) is, 
using (8.2), equal to Sh(t), and we therefore have, that 

S G (t)=S G (ntWt) + S h (t) (8.9) 

Eq. (8.9), which is valid for all macroscopic systems in LTE, gives an important connec- 
tion between Sq and Sh for such systems. 

For an isolated system with a given Hamiltonian H, we write X = Th{X) for the flow 
in the phase space. The ensemble density /i(X,t) then evolves according to the Liouville 
equation, 

= £ HfI = -div(^F fl ) = -V/x • F H , (8.10) 

subject to some initial condition fi(x,0) = pio(X). (We think of X, J 7 , V as (2dN)- 
dimensional vectors in the phase space and in the last equality have used the fact that 
divJ-H = 0.) As was already known to Gibbs, and is proven in almost every text book 
on statistical mechanics, Sa(t) is constant in time for an isolated system. Hence, using 
(8.9) 

s G (t) = s G (m\v t ) + s h (t) = o (8.ii) 

or 

SaifhWt) = ~S h (t) = -R^t) (8.12) 

Note that in going from (8.11) to (8.12) we are glossing over the difference between 
microscopic and macroscopic time scales, see [11, 25, 34]. In the same spirit we will also 
have that Ri(t) = (Ri) where Ri is a suitable microscopic function and the average is 
with respect to fi t , not v t . 

As t — > oo, we expect that fi t and u t will both "approach" the Gibbs distribution 
/i cq appropriate for a macroscopic system in equilibrium with a uniform temperature 
and density. Note however, that while Shit) approaches S cq , the Gibbs entropy, Sa(ii t ), 
being constant, clearly does not. Hence if fJ-o(X) is not an equilibrium state, the limit, as 
t — > oo, of Sci/J'tWt) will be Sg(0) — S eq , which is negative and can be large, i.e. extensive, 
indicating the "entanglement" of ^t(X) necessary to keep So constant. 

An explicit example where Sc(^tWt) can be shown to be extensive for large t, occurs 
for non interacting point particles moving among a periodic set of fixed scatterers, the 
periodic Lorentz gas or Sinai billiard, which are started at t = in a product measure 
with the same energy (speed) and a macroscopically nonuniform density. The density 
then evolves on the macroscopic scale according to the diffusion equation, see [11, 25], 
and references there. This leads to a uniformization of the density and an increase in 
Sh{t) while Sa(t) remains constant. The entanglement of ^t(X) corresponds here to 
the build up of correlations between the positions and velocities of the particles in the 
system as it evolves towards a spatially uniform state. (The lack of interaction between 
the particles makes this system a bit unphysical.) 
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We turn now to the statistical mechanics of open systems, e.g. a fluid with shear 
flow, whose macroscopic behavior is described by the compressible Navier-Stokes eqs. 
Modeling such a system microscopically as one driven by stochastic thermal reservoirs 
and some external forces, e.g. rough walls with constant temperatures and velocities, as 
in section 4b, the time evolution of n(X, t) will now be given by [4,5] 

^ = (C H + C cx + J2£«)v (8-13) 

Here K, a represents the stochastic effect of the ath reservoir which tries to bring the 
system to equilibrium with a temperature T a , and £ cx fi = — jF ex VyU, assuming that the 
external force jF ex is (phase-space) divergence free. The Gibbs entropy Sa(t) will no 
longer be constant in time so (8.12) will no longer hold. Assuming the system to be in 
LTE, we shall have instead of (8.12) the behavior given in (8.6), with 

S h (t) = (R l )-Y,(Ja)/T a (8.14) 

where (•) is the average with respect to fj,(X,t) of phase space functions Ri and J a ; J a 
specifies the rate at which energy flows from the system to the ath reservoir. As in 
(8.6) we are assuming here that in our system, the spatial region in contact with the ath 
reservoir is itself at a temperature very close to T a . (Compare (8.14) with (4.3) where 
there is no assumption of LTE; 1Z there corresponds to Sa{nt\vt) + Ri)- 

We may assume that our ensemble density /i(X,t), evolving according to (8.13) will, 
in contrast to what happens in an isolated system, approach a smooth stationary density 
fl(X) while v t {X) -> V(X) = u(X\M p ) [4,5,10]. Thus for t -> oo, S h (t) will approach 
the hydrodynamic entropy Sh in the MSNS, which according to (8.2) is given by 

S h = - J u(X)\ogu(X)dX. (8.15) 

Setting Ss(t) = 0, in (8.14), and letting Ri be the internal entropy production in the 
stationary state, we will have 

Ri = J2( J «)fi/T a - (8.16) 

Comparing (8.16) with (8.7) we can identify the average of J a with respect to fx, {J a )p., 
with the hydrodynamic steady state energy flux J a , which, for a system in LTE, should 
be just (J a )p- Since we further expect that the stochasticity of the reservoirs will make 
the stationary fx a "smooth" density, independent of the initial fj, , we should also have 
here 

Sc(t) —> 0, and Sai^tl^i) ~^ Sg(Jx\v) for t — > oo (8-17) 

Let us turn now to the microscopic modeling of the open macroscopic system by "ther- 
mostatted" deterministic forces, J- ts (X), which conserve the energy but whose divergence 
does not vanish, e.g. the Maxwell demon boundary drives discussed in this paper. For 
the sake of simplicity we shall treat !Fx, s (X) as if it was smooth, but think of it as acting 
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only in the vicinity of the boundaries; the transition to a boundary term should then be 

possible, see Appendix. The time evolution of a microstate X will now be given by the 
deterministic equation 

X = T H {X)+F*{X), (8.18) 

and that of fi(X, t) by 

^ = -div[(^ + ^ts)/i]. (8.19) 
Taking the time derivative of Sa{t), we get, as in (4.6), 

S G (t) = J m(X) divT ts (X) dX = -M(t) (8.20) 

We also have, using (8.8) and (8.9), for systems in LTE, that 

S h (t) = j J Ht(X) log v t (X) dx = - J F H {V ■ fi t ) log v t dx 

- J div(T ts )fi t logv t dx- J fi t (X)^- t \ogu t (X)dX (8.21) 

In the last term on the rhs of (8.21) we can replace fi t by v t since the time derivative of 
log v t at a fixed X, involves only phase space functions whose expectation corresponds to 
the hydrodynamic variables. After this replacement the last term becomes ^ / v t (X) dX, 
which vanishes by normalization. We are thus left with 

S h (t)= - j{F H -Vfi t )\ogv t dX 

+ j ' titFto-VlogvtdX (8.22) 

We now want to argue that for "smooth" thermostatted forces we can again replace 
fit by v t in the last term in (8.22), obtaining 

J v t F ts ■ V log v t dX = - j F ts ■ Vv t dX 

= Ju t div^ ts dX ee -Mt (t) (8.23) 

We also argue that the first term on the right side of (8.22) which is the only term 
present in an isolated system, is just Ri(t). Accepting these "arguments" we finally get 
for systems with thermostatted forces in LTE, 

S h (t) = R l (t)-M l (t), (8.24) 

We note that (8.24) gives a decomposition of the rate of change of Sh into an internal 
part, R[, coming from the dissipative fluxes inside the system, and a thermostatted part, 
— Mi, coming from the thermostat forces jF ts . Comparing this with (8.4) and remem- 
bering that J-'ts does not change n or e, M/ will contribute only to the last term there, 
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— / ^^- t (\nu 2 )dr. In fact, Mi corresponds to the surface term in (8.5), which for the 
models considered in this paper is equivalent to the second term on the right side of 
(8.6): it represents the production (or reduction) of Sh due to the conversion of thermal 
into directed energy by the Maxwell demons producing the jF ts at the boundary. 

Interestingly, the contribution of the Tt s to Sh, —Mi(t), is just the rate of change of 
Sa(t) when ji t = u t . It is tempting to try and bypass the formal manipulations leading 
to (8.24) and derive that relation more directly from the definition of Sb{M) in (8.3) as 
logT(A / f), but we have not succeeded in doing this so far. 

Waiting now for a time t which is long enough for the system to become approximately 
stationary on the macroscopic level, we would have Shit) ~ 0, for t > t , and thus 

Ri = M t (8.25) 

Eq. (8.25) explains the equality between (5.12) and (5.15). If it is further true that divX 
is sufficiently smooth for its average to be well approximated by Mi, then we would have 
M — Mi — Ri. This appears to be the case in many situations, including the Maxwell 
demon boundaries considered in this paper, where divX is an additive function of the 
coordinates and velocities of each particle so that its average, M, depends only on the 
one particle distribution function at the wall. This leads to the equality, M Ri, which 
we observe in our simulations. We actually expect the equality (8.25) to hold for general 
thermostatted SNS as long as LTE holds in the SNS. It should in particular hold for 
macroscopic fluids, driven by rules (2.3) or (2.4) even when the flow is no longer laminar. 
We hope to test this via simulations on large systems. 
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Appendix. Special reflection rules 

It seems reasonable to expect that time reversible reflection rules, like rule b, can be 
obtained as limits of the usual smooth Gaussian thermostatted dynamics [7, 8]. We 
describe here a very simple example of such a limit which also provides an illustration of 
the relation between Mi and R t , discussed in section 8. 

Let a constant oblique force F = (E cos (3, — E sin/3) act in the half-plane y > L/2, 
above the box, and a symmetric force F = (— E cos/5, Esin/3) act below the box, in the 
area y < —L/2. Here < (3 < 7r/2isa fixed parameter, and E > is very large, 
'almost infinite'. Such a force effectively replaces the walls. The force is accompanied 
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by a Gaussian constraint which keeps the kinetic energy of the particle fixed [7, 8]. The 
motion of a single particle in the region y > L/2 is then governed by the equations, 

x = v x , y = v y 

and v = Tts-, is given by 



v x 



where 



a: 



= E cos ft — av x 
= — E sin ft — av y , 

Ev x cos ft — Ev y sin ft 
v 2 + v 2 

x ' y 



(A.l) 



so that v 2 is constant. Symmetric equations hold for y < —L/2. 

Taking now the limit E — > oo this model reduces to our Maxwell demon model (2.1) 
with a specific function f = fp. To obtain that, we take advantage of the conservation 
of kinetic energy and rewrite the system (A.l) as 



x 

y 

6 



v cos 6 
v sin 9 

--sin(e + ft) 



(A.2) 



where 9 = ta,n~ 1 (v y /v x ) is the angle between the velocity vector (v x ,v y ) and the x axis. 
The last equation in ( |A.2| ) is independent of the first two, and has an implicit solution 
given by 



In 



1 + cos(# + ft) 



2™ 

V 



(A.3) 



l-cos(9 + ft)_ 

It is, however, more convenient to differentiate 9 with respect to the height, h = y — L/2, 
yielding 

E sin(fl + ft) 



d9 _ 

dh v 

An implicit solution of this equation is 



sin i 



(A.4) 



Eh 

(9 + ft)-cosft-smft-\n\sin{9 + ft)\ = T 



(A.5) 



Since h = both as the particle enters the force zone and as it leaves it, the relation 
between the incoming angle, = 9 in and the outgoing angle, ip = —9 out , is determined 
by the equation 

F p {-i>) = F p {ip) (A.6) 

where 

F (9) := (9 + ft)- cos ft - sin ft ■ In | sin(fl + ft) | (A.7) 
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These equations do not contain the force strength, E. By taking the limit E — > oo 
we simply assure that the particle leaves the force zone instantly, the moment it enters 
it. Therefore, the infinite force acts like a wall at which the particle gets reflected, 
and the outgoing angle ip is related to the incoming angle <p by equation (2.1) with 
f fi (<p) = -Ff 1 F fl (<p). 

For E very large the time any particle spends in the force zone is extremely short, 
so we can neglect possible collisions between particles while one of them is in that zone. 
In fact if a particle enters the force zone at time s and leaves it at time s + r we can 
assume that there is no other particle in the force zone between s and s + r. Under these 
conditions the analysis of section 8 takes on a particularly simple form. 

Assuming that the stationary state is one of LTE, the one particle distribution in the 
region of the field is a local Maxwellian where \P(r, v) 

tt(r,v) = n(r)(27rT„,)- 1 exp{-(v- v w ) 2 /2T w } (A.8) 

Using the definition of M; in (8.23) with jF ts , given in (A.l), we obtain 

Mi = -Jdv j(-^-F ts )^{v,v)dv 

= ~ [dr [ { ^M-^(v,v)dv 

J J 1 w 

= ^ / dr f v^(r,v)dv (A.9) 

In going from the second to the third equality on the right side of ( |A.9|) we have used 
the conservation of energy by ^tsj i-e. v ■ JF ts = 0. The final term is easily recognized 
as corresponding to the last term in (8.4), giving the entropy production due to JF ts . It 
becomes equal to Ri in (4.2) in the limit E —>■ oo, giving Mi = Ri, for such systems in 
LTE. 
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